MODULE SURFWS_INIT_MLOFF_MOD
CONTAINS

SUBROUTINE SURFWS_INIT_MLOFF(KIDIA, KFDIA, KLON, KLEVSN, NCL, PMU0,PSDOR,        &  ! Input
                     &  PTSOIL, PTSKIN, &
                     &  ZDSNTOT, ZSNDEPTH,                  &       
                     &  ZSNPERT,                                           &  ! Input
                     &  ZDSNR,PTSN, PRSN, PSSN, PWSN,PASN,                 &  ! Input
                     &  PTSNWS,PSSNWS,PRSNWS,PWSNWS,                       & ! Output
                     &  PTSNTOP, PTSNBOTTOM, PTSNMIDDLE,                   & ! Output
                     &  PRSNMAX, ZSADEPTH, KLEVSNA, KLEVMID, ZACTDEPTH,    & ! Output
                     &  PTMINCL,PRMINCL,                                   & ! Output
                     &  PTCONSTAVG, PTCONSTSTD,                            & ! Output
                     &  PRCONSTAVG, PRCONSTSTD, PRSNTOP,                   & ! Output
                     &  YDCST, YDSOIL )

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST  , ONLY : TCST
USE YOS_SOIL , ONLY : TSOIL

USE ABORT_SURF_MOD

! (C) Copyright 2017- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**** *SURFWS_INIT* - Snow warm start multi-layer 
!     PURPOSE.
!     --------
!          THIS ROUTINE CONTROLS THE WARM START OF MULTI-LAYER SCHEME IN
!          FC MODE

!**   INTERFACE.
!     ----------
!          *SURFWS_CTL* IS CALLED FROM *SRFSN_DRIVER*.

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!    *KIDIA*      START POINT


!     INPUT PARAMETERS (REAL):
!    *PTMST*      TIME STEP                                      S

!     INPUT PARAMETERS (LOGICAL):
!    *LDLAND*     LAND/SEA MASK (TRUE/FALSE)

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!    *PWLM1M*     SKIN RESERVOIR WATER CONTENT                   kg/m**2


!     OUTPUT PARAMETERS AT T+1 (UNFILTERED,REAL):
!    *PWL*        SKIN RESERVOIR WATER CONTENT                   kg/m**2


!     OUTPUT PARAMETERS (DIAGNOSTIC):

!    *PDHIIS*     Diagnostic array for interception layer (see module yomcdh)

!     METHOD.
!     -------
!          

!     EXTERNALS.
!     ----------
!          NONE.

!     REFERENCE.
!     ----------
!          

!     Modifications:
!     Original   G. Arduini      ECMWF     28/07/2017

!     ------------------------------------------------------------------

IMPLICIT NONE

! Declaration of arguments 
INTEGER(KIND=JPIM), INTENT(IN) :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN) :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN) :: KLON
INTEGER(KIND=JPIM), INTENT(IN) :: KLEVSN
INTEGER(KIND=JPIM), INTENT(IN) :: NCL

REAL(KIND=JPRB), INTENT(IN)    :: ZDSNTOT(:), PTSN(:), PRSN(:), PSSN(:), PWSN(:)
REAL(KIND=JPRB), INTENT(IN)    :: ZSNDEPTH(:,:)
REAL(KIND=JPRB), INTENT(IN)    :: ZDSNR(:,:)
REAL(KIND=JPRB), INTENT(IN)    :: PSDOR(:)

REAL(KIND=JPRB), INTENT(IN)    :: PMU0(:)
REAL(KIND=JPRB), INTENT(IN)    :: PTSOIL(:)
REAL(KIND=JPRB), INTENT(IN)    :: PTSKIN(:)
REAL(KIND=JPRB), INTENT(IN)    :: PASN(:)

TYPE(TCST)     , INTENT(IN) :: YDCST
TYPE(TSOIL)    , INTENT(IN) :: YDSOIL

! Output Variables
REAL(KIND=JPRB), INTENT(OUT)    :: PTSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PSSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PRSNWS(:,:)
REAL(KIND=JPRB), INTENT(OUT)    :: PWSNWS(:,:)

REAL(KIND=JPRB), INTENT(OUT)    :: PTSNTOP(:), PTSNBOTTOM(:), PTSNMIDDLE(:)
REAL(KIND=JPRB), INTENT(OUT)    :: PRSNTOP(:)

REAL(KIND=JPRB), INTENT(OUT)    :: PRSNMAX(:)
REAL(KIND=JPRB), INTENT(OUT)    :: ZSADEPTH(:)

! Active number of layers:
INTEGER(KIND=JPIM),INTENT(INOUT)  :: KLEVSNA(:)

! Active depth layer (glaciers):
REAL(KIND=JPRB)    :: ZACTDEPTH(:)

! Cluster values:
INTEGER(KIND=JPIM),            INTENT(OUT)  :: KLEVMID(:)
INTEGER(KIND=JPIM),            INTENT(OUT)  :: PTMINCL(:), PRMINCL(:) 
REAL(KIND=JPRB),               INTENT(OUT)  :: PRCONSTAVG(:,:), PRCONSTSTD(:,:)
REAL(KIND=JPRB),               INTENT(OUT)  :: PTCONSTAVG(:,:), PTCONSTSTD(:,:)

REAL(KIND=JPRB)                 :: ZSNPERT

! Local variables:
INTEGER(KIND=JPIM)              :: JL,JK,KNACC
INTEGER(KIND=JPIM)              :: KCL
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTDIST, ZRDIST
REAL(KIND=JPRB)                 :: ZFEATA, ZFEATB, ZFEATC 

! Look-up tables for warm-start
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCENTRADAY2, ZTCENTRBDAY2 ,ZTCENTRCDAY2 ,ZTCENTRADAY3 ,ZTCENTRBDAY3 ,ZTCENTRCDAY3 ,&
                                    ZTCENTRADAY4 ,ZTCENTRBDAY4 ,ZTCENTRCDAY4 ,ZTCENTRADAY5 ,ZTCENTRBDAY5 ,ZTCENTRCDAY5 ,&
                                    ZTCENTRADAY5M ,ZTCENTRBDAY5M ,ZTCENTRCDAY5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCENTRANIGHT2 ,ZTCENTRBNIGHT2 ,ZTCENTRCNIGHT2 ,ZTCENTRANIGHT3 ,&
                                    ZTCENTRBNIGHT3 ,ZTCENTRCNIGHT3 ,ZTCENTRANIGHT4 ,ZTCENTRBNIGHT4 ,&
                                    ZTCENTRCNIGHT4 ,ZTCENTRANIGHT5 ,ZTCENTRBNIGHT5 ,ZTCENTRCNIGHT5 ,&
                                    ZTCENTRANIGHT5M ,ZTCENTRBNIGHT5M ,ZTCENTRCNIGHT5M 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVGDAY2 ,ZTCONSTSTDDAY2 ,ZTCONSTAVGDAY3 ,ZTCONSTSTDDAY3 ,&
                                    ZTCONSTAVGDAY4 ,ZTCONSTSTDDAY4 ,ZTCONSTAVGDAY5 ,ZTCONSTSTDDAY5 ,&
                                    ZTCONSTAVGDAY5M ,ZTCONSTSTDDAY5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVGNIGHT2 ,ZTCONSTSTDNIGHT2 ,ZTCONSTAVGNIGHT3 ,ZTCONSTSTDNIGHT3,&
                                    ZTCONSTAVGNIGHT4 , ZTCONSTSTDNIGHT4 ,ZTCONSTAVGNIGHT5 ,ZTCONSTSTDNIGHT5,&
                                    ZTCONSTAVGNIGHT5M ,ZTCONSTSTDNIGHT5M

REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTCONSTAVG, ZRCONSTAVG, ZTCONSTSTD, ZRCONSTSTD
REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZTCENTRA, ZTCENTRB, ZTCENTRC
REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZRCENTRA, ZRCENTRB, ZRCENTRC


REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRIDAY2 ,ZTMLRADAY2 ,ZTMLRBDAY2 ,ZTMLRCDAY2 ,ZTMLRDDAY2 ,ZTMLREDAY2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRINIGHT2 ,ZTMLRANIGHT2 ,ZTMLRBNIGHT2 ,ZTMLRCNIGHT2 ,ZTMLRDNIGHT2 ,ZTMLRENIGHT2
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRIDAY3 ,ZTMLRADAY3 ,ZTMLRBDAY3 ,ZTMLRCDAY3 ,ZTMLRDDAY3 ,ZTMLREDAY3
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZTMLRINIGHT3 ,ZTMLRANIGHT3 ,ZTMLRBNIGHT3 ,ZTMLRCNIGHT3 ,ZTMLRDNIGHT3 ,ZTMLRENIGHT3
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA2 ,ZRCENTRB2 ,ZRCENTRC2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA3 ,ZRCENTRB3 ,ZRCENTRC3 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA4 ,ZRCENTRB4 ,ZRCENTRC4
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA5 ,ZRCENTRB5 ,ZRCENTRC5
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCENTRA5M ,ZRCENTRB5M ,ZRCENTRC5M
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRCONSTAVG2 ,ZRCONSTSTD2 ,ZRCONSTAVG3 ,ZRCONSTSTD3 ,&
                                    ZRCONSTAVG4 ,ZRCONSTSTD4 ,ZRCONSTAVG5 ,ZRCONSTSTD5 ,ZRCONSTAVG5M ,ZRCONSTSTD5M 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI2 ,ZRMLRA2 ,ZRMLRB2 ,ZRMLRC2 ,ZRMLRD2 ,ZRMLRE2 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI3 ,ZRMLRA3 ,ZRMLRB3 ,ZRMLRC3 ,ZRMLRD3 ,ZRMLRE3 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI4 ,ZRMLRA4 ,ZRMLRB4 ,ZRMLRC4 ,ZRMLRD4 ,ZRMLRE4 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI5 ,ZRMLRA5 ,ZRMLRB5 ,ZRMLRC5 ,ZRMLRD5 ,ZRMLRE5 
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI5M ,ZRMLRA5M ,ZRMLRB5M ,ZRMLRC5M ,ZRMLRD5M ,ZRMLRE5M

REAL(KIND=JPRB), DIMENSION(NCL)  ::  ZTMLRA, ZTMLRB, ZTMLRC, ZTMLRD, ZTMLRE, ZTMLRI
REAL(KIND=JPRB), DIMENSION(NCL)  :: ZRMLRI, ZRMLRA, ZRMLRB, ZRMLRC, ZRMLRD, ZRMLRE
REAL(KIND=JPRB)                  :: ZP1, ZP2, ZP3, ZP4, ZP5
REAL(KIND=JPRB)                  :: ZSDORTHR
REAL(KIND=JPRB)                  :: ZEPSILON

REAL(KIND=JPHOOK)                  :: ZHOOK_HANDLE


!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFWS_INIT_MLOFF_MOD:SURFWS_INIT_MLOFF',0,ZHOOK_HANDLE)

!    -----------------------------------------------------------------

ASSOCIATE(RTT=>YDCST%RTT,RLMLT=>YDCST%RLMLT, RPI=>YDCST%RPI,    &
        & RLWCSWEA=>YDSOIL%RLWCSWEA, RLWCSWEB=>YDSOIL%RLWCSWEB, &
        & RLWCSWEC=>YDSOIL%RLWCSWEC, RTEMPAMP=>YDSOIL%RTEMPAMP, &
        & RDSNMAX=>YDSOIL%RDSNMAX, RHOMINSND=>YDSOIL%RHOMINSND, &
        & RHOMAXSN_NEW=>YDSOIL%RHOMAXSN_NEW, RDAT=>YDSOIL%RDAT  )

! mid-point first soil level
ZSADEPTH(:)=0.5_JPRB*RDAT(1)
ZSDORTHR=50._JPRB

ZEPSILON  = 10E4*EPSILON(ZEPSILON)

! 0.1 Define centroids 
! A: soT-skt B: soT-snT ; C: rsn/rsnmin
ZTCENTRADAY2(1:NCL)=(/ 10.665,-0.746,-11.857,7.005,20.408,-4.571,-0.414,1.902,12.724,7.379,24.801,-9.194,-15.782,2.555,5.719,-5.452,3.056,17.634,-7.019,-2.571,-2.65,3.952,15.21,0.762,8.964,9.363,4.337  /)
    ZTCENTRBDAY2(1:NCL)=(/ 4.781,-0.267,-1.156,2.776,8.314,-1.646,0.063,2.046,5.265,3.294,9.803,-1.012,-1.084,0.515,2.159,-0.613,1.134,7.235,-1.633,-0.1,-0.92,1.695,6.044,0.423,9.546,3.372,3.996  /)
    ZTCENTRCDAY2(1:NCL)=(/ 2.494,2.595,4.17,5.094,2.082,2.899,5.544,2.275,2.339,2.231,1.962,4.18,4.105,2.342,2.275,4.885,4.896,2.155,3.299,5.529,2.843,2.294,2.205,2.509,1.935,2.248,2.103  /)
    ZTCENTRADAY3(1:NCL)=(/ 3.681,18.971,-7.308,11.437,0.007,16.422,24.398,6.981,-0.621,-13.271,13.776,-4.043,5.222,-9.71,33.778,9.633,2.111,-1.882,1.674,28.233,-6.001,2.332,7.409,-2.863,16.404,21.543,9.165  /)
    ZTCENTRBDAY3(1:NCL)=(/ 1.192,7.303,-1.917,4.373,-0.11,6.485,9.431,2.607,-0.017,-1.34,5.264,-1.859,1.918,-1.341,12.739,2.014,2.218,-1.048,0.498,10.614,-1.13,0.585,2.138,-0.295,4.092,8.052,3.804  /)
    ZTCENTRCDAY3(1:NCL)=(/ 2.402,2.155,3.088,2.423,2.713,2.235,2.124,2.307,5.504,4.062,2.335,2.788,2.35,3.974,1.913,2.228,2.292,2.84,2.427,2.007,4.273,4.91,5.157,5.258,1.887,2.144,2.445  /)
    ZTCENTRADAY4(1:NCL)=(/ 10.56,0.149,19.497,-9.183,3.922,35.092,16.334,3.376,-1.685,25.12,12.907,1.618,8.435,-1.946,6.902,17.105,-3.95,9.611,14.728,12.546,5.064,28.846,19.603,-6.319,-0.144,-12.722,22.192  /)
    ZTCENTRBDAY4(1:NCL)=(/ 3.897,-0.062,7.189,-1.702,0.519,11.943,3.312,1.252,-0.882,8.776,2.498,0.892,3.351,-0.026,2.114,6.206,-1.564,1.795,5.447,4.808,1.871,9.923,4.285,-1.773,0.064,-1.04,7.703  /)
    ZTCENTRCDAY4(1:NCL)=(/ 2.596,2.81,2.3,3.588,5.09,2.175,2.107,2.485,2.878,2.26,2.292,2.541,2.568,5.43,2.593,2.339,3.21,2.714,2.489,2.522,2.576,2.273,1.983,3.4,5.285,4.341,2.261  /)
    ZTCENTRADAY5(1:NCL)=(/ 32.408,6.341,-4.106,18.849,11.891,-11.135,-0.284,0.342,28.252,24.528,9.695,-2.07,7.203,-7.373,39.383,-1.44,3.577,14.854,3.385,5.561,7.991,14.081,21.734,10.979,18.98,1.741,16.374  /)
    ZTCENTRBDAY5(1:NCL)=(/ 8.96,3.454,-1.153,6.721,4.02,-1.178,0.145,0.195,7.992,7.344,3.599,0.2,0.071,-1.379,10.499,-0.702,1.948,2.654,0.051,1.695,2.396,4.975,6.354,1.03,4.506,1.251,5.356  /)
    ZTCENTRCDAY5(1:NCL)=(/ 2.403,2.48,3.326,2.347,2.571,3.956,5.006,2.974,2.425,2.408,2.58,5.37,4.161,3.54,2.379,3.426,2.587,2.853,3.927,2.73,2.732,2.477,2.421,3.58,2.545,2.665,2.448  /)
    ZTCENTRADAY5M(1:NCL)=(/ 4.785,19.234,0.792,8.587,26.497,13.488,6.71,12.668,-0.841,37.772,16.331,22.685,2.224,0.08,31.056,3.944,12.227,24.207,3.064,1.24,15.649,19.373,9.904,-1.57,7.331,10.082,5.897  /)
    ZTCENTRBDAY5M(1:NCL)=(/ 0.591,1.374,0.773,2.902,5.902,5.384,4.648,0.698,-0.028,7.553,4.797,5.595,0.503,0.49,6.097,3.673,3.204,1.952,1.285,2.971,1.096,5.13,0.77,0.411,0.659,5.089,2.142  /)
    ZTCENTRCDAY5M(1:NCL)=(/ 4.257,3.848,3.403,3.139,3.028,2.929,3.068,4.158,3.774,2.891,3.089,3.018,4.565,4.746,3.064,3.206,3.231,3.753,3.189,3.483,3.932,3.062,4.097,5.224,4.193,2.922,3.071  /)


    ZTCENTRANIGHT2(1:NCL)=(/ -0.651,14.093,6.256,9.776,19.373,-8.883,3.688,26.036,0.561,7.497,-5.54,17.692,0.838,5.362,22.849,3.867,10.815,30.659,15.924,3.715,21.54,-2.544,12.537,8.269,-2.698,2.156,7.078  /)
    ZTCENTRBNIGHT2(1:NCL)=(/ -0.186,6.389,3.488,4.084,8.645,-1.0,1.561,10.975,0.084,2.799,-1.574,7.641,0.511,2.239,9.833,3.08,5.212,11.767,6.944,0.616,7.849,-0.782,5.435,4.158,-0.997,1.358,2.534  /)
    ZTCENTRCNIGHT2(1:NCL)=(/ 2.432,2.2,2.179,2.264,2.111,4.34,2.241,1.945,5.193,2.196,3.266,2.151,2.307,2.302,2.009,2.201,2.388,1.949,2.178,5.111,1.856,4.485,2.223,2.336,2.441,2.373,4.791  /)
    ZTCENTRANIGHT3(1:NCL)=(/ 6.557,21.885,-0.557,12.786,30.201,16.116,9.468,2.275,17.978,-3.906,24.217,7.965,0.38,10.992,3.268,-8.062,-4.655,14.491,35.26,19.872,3.779,0.939,5.143,9.434,26.98,-2.078,6.393  /)
    ZTCENTRBNIGHT3(1:NCL)=(/ 2.757,8.877,-0.337,5.145,11.772,6.608,2.722,1.255,7.294,-0.915,9.685,3.386,-0.055,4.611,0.156,-1.426,-1.748,5.851,13.022,8.135,1.688,0.524,2.325,3.996,10.567,-1.156,0.624  /)
    ZTCENTRCNIGHT3(1:NCL)=(/ 2.32,2.118,2.533,2.349,2.025,2.172,4.207,2.368,2.171,4.513,2.06,2.294,5.187,2.278,5.071,3.454,2.597,2.18,1.883,2.179,2.31,2.495,2.267,2.25,1.976,3.035,4.949  /)
    ZTCENTRANIGHT4(1:NCL)=(/ 10.346,18.492,1.675,30.468,-3.631,11.878,22.46,7.479,42.471,16.697,-7.42,0.183,27.332,7.652,0.331,13.44,3.599,34.273,6.228,-1.455,4.756,11.882,9.075,20.353,15.044,24.729,3.274  /)
    ZTCENTRBNIGHT4(1:NCL)=(/ 3.969,6.637,0.823,10.492,-1.499,4.469,7.888,0.521,14.675,6.084,-1.503,-0.123,9.37,2.939,0.001,4.912,0.012,11.747,2.403,-1.125,1.895,2.52,3.261,7.353,5.421,8.676,1.346  /)
    ZTCENTRCNIGHT4(1:NCL)=(/ 2.451,2.285,2.456,2.129,3.237,2.427,2.236,4.481,1.999,2.276,3.316,2.533,2.169,2.463,5.218,2.43,4.848,2.106,2.452,2.669,2.469,3.239,2.6,2.239,2.418,2.235,2.42  /)
    ZTCENTRANIGHT5(1:NCL)=(/ 3.062,20.714,12.436,30.858,1.084,21.748,17.985,6.939,27.569,14.264,5.76,11.602,-5.275,10.587,8.819,3.14,34.898,14.982,24.912,8.65,0.465,4.99,41.511,16.229,18.683,-1.105,23.734  /)
    ZTCENTRBNIGHT5(1:NCL)=(/ 1.318,5.38,4.062,8.381,0.49,7.329,4.37,2.356,7.392,4.859,0.114,1.183,-1.27,3.606,2.901,-0.189,9.466,2.604,8.142,0.324,0.134,1.924,10.709,5.419,6.479,-0.929,5.933  /)
    ZTCENTRCNIGHT5(1:NCL)=(/ 2.542,2.471,2.516,2.336,2.672,2.287,2.674,2.561,2.394,2.424,4.304,3.558,3.697,2.497,2.583,4.339,2.348,3.061,2.337,4.267,4.969,2.505,2.301,2.406,2.316,3.132,2.434  /)
    ZTCENTRANIGHT5M(1:NCL)=(/ 8.979,19.759,0.575,29.292,15.061,14.129,6.453,40.803,15.995,10.028,23.28,2.944,3.806,25.493,10.616,7.739,5.169,12.397,33.645,1.44,6.382,21.858,26.95,18.062,21.286,12.538,17.699  /)
    ZTCENTRBNIGHT5M(1:NCL)=(/ 2.895,4.031,0.39,6.091,0.847,5.357,1.434,8.001,3.665,0.703,4.381,0.51,2.021,6.223,4.952,0.553,0.531,0.712,6.638,1.028,3.973,1.329,2.213,1.076,6.15,3.106,5.587  /)
    ZTCENTRCNIGHT5M(1:NCL)=(/ 3.097,3.375,4.654,3.04,4.053,2.851,3.108,2.855,3.303,4.199,3.333,4.552,2.971,2.913,2.867,4.519,4.497,4.146,2.957,3.127,2.969,3.854,3.697,3.948,2.822,3.271,2.905  /)


    ZTCONSTAVGDAY2(1:NCL)=(/ -3.582,-0.606,-9.657,-3.525,-3.362,-7.4,-0.184,-0.975,-3.877,-3.46,-3.249,-2.461,-10.347,-0.349,-2.469,-0.23,-1.44,-3.193,-8.025,-0.436,-1.448,-1.786,-3.196,-0.416,-1.276,-3.102,-2.183  /)
    ZTCONSTSTDDAY2(1:NCL)=(/ 5.232,6.754,43.162,6.192,4.638,31.216,0.743,2.813,4.863,5.17,4.409,18.944,40.037,2.493,4.609,1.401,4.468,4.479,34.157,5.73,11.595,3.952,4.368,3.197,2.775,5.037,4.191  /)
    ZTCONSTAVGDAY3(1:NCL)=(/ -1.165,-2.648,-15.721,-2.772,-0.062,-2.982,-2.471,-2.355,0.037,-12.187,-2.823,-3.115,-1.985,-13.782,-2.717,-1.607,-0.602,-0.018,-0.316,-2.034,-0.958,-0.415,-2.728,-0.112,-1.877,-2.739,-2.726  /)
    ZTCONSTSTDDAY3(1:NCL)=(/ 5.033,3.306,53.491,3.839,3.903,4.427,3.166,5.966,3.275,45.145,3.914,21.061,5.299,50.76,3.023,5.268,1.86,3.696,2.995,2.906,9.397,3.514,4.51,1.029,3.428,3.317,3.962  /)
    ZTCONSTAVGDAY4(1:NCL)=(/ 4.81,5.346,4.535,8.548,8.029,4.787,12.482,5.887,14.009,5.53,11.806,4.466,5.635,10.959,8.123,4.677,28.795,11.503,5.188,4.517,6.849,5.581,8.354,15.499,2.262,10.025,5.326  /)
    ZTCONSTSTDDAY4(1:NCL)=(/ 7.59,34.417,4.353,6.312,19.905,4.314,32.522,29.629,60.524,4.585,32.619,29.098,23.074,41.903,33.539,4.448,80.154,30.674,4.682,4.417,33.658,4.317,9.189,46.968,18.749,5.672,4.467  /)
    ZTCONSTAVGDAY5(1:NCL)=(/ 5.547,2.93,20.857,4.186,5.055,5.396,0.867,5.216,5.423,5.32,3.76,3.649,10.37,5.04,5.843,7.831,2.108,8.869,9.39,6.521,6.846,4.814,5.505,9.86,7.911,3.371,5.243  /)
    ZTCONSTSTDDAY5(1:NCL)=(/ 2.837,24.369,69.91,2.689,13.499,3.59,9.121,32.252,2.744,2.955,9.7,18.084,16.791,4.261,2.377,43.162,17.515,7.799,30.686,29.127,26.947,15.207,2.998,6.299,17.192,25.22,9.313  /)
    ZTCONSTAVGDAY5M(1:NCL)=(/ 5.646,6.23,2.056,1.511,4.194,0.762,2.103,7.308,1.35,4.683,2.652,3.742,5.821,0.752,4.599,12.78,2.731,5.607,1.733,24.743,6.693,3.278,7.089,1.961,6.919,-0.456,0.945  /)
    ZTCONSTSTDDAY5M(1:NCL)=(/ 14.428,2.949,15.951,9.658,1.265,1.331,22.858,6.465,10.196,1.146,1.72,1.456,18.831,4.874,1.099,46.176,2.18,1.448,12.549,56.278,4.953,1.647,9.219,1.756,11.429,1.163,7.715  /)
    ZTCONSTAVGNIGHT2(1:NCL)=(/ -0.671,-3.784,-3.813,-3.622,-4.092,-14.571,-2.602,-3.072,-0.222,-3.638,-8.44,-3.791,-0.522,-3.319,-3.16,-2.566,-4.692,-2.332,-3.889,-1.004,-2.729,-1.155,-3.496,-3.708,-2.295,-1.425,-3.144  /)
    ZTCONSTSTDNIGHT2(1:NCL)=(/ 8.819,5.141,5.574,5.019,5.324,56.195,5.28,4.579,2.074,5.493,35.837,4.907,2.572,5.369,4.614,4.895,5.82,3.46,5.011,3.329,3.559,11.968,4.829,5.359,14.389,4.055,5.614  /)
    ZTCONSTAVGNIGHT3(1:NCL)=(/ -2.434,-2.655,-0.446,-2.693,-1.988,-2.758,-2.188,-0.925,-2.763,-0.353,-2.543,-2.479,-0.146,-2.826,-0.387,-13.81,-4.236,-3.094,-1.901,-2.972,-1.673,-0.333,-2.277,-2.751,-2.329,-2.66,-0.829  /)
    ZTCONSTSTDNIGHT3(1:NCL)=(/ 4.066,4.243,7.907,3.453,2.817,3.657,3.176,3.154,3.267,1.793,3.098,4.037,2.312,3.868,1.682,55.968,27.478,5.146,2.688,3.83,3.904,1.771,4.412,3.694,3.089,21.875,1.938  /)
    ZTCONSTAVGNIGHT4(1:NCL)=(/ 3.812,4.609,2.846,5.016,23.403,4.233,4.851,9.914,4.769,4.96,8.52,4.016,4.701,3.538,2.031,4.573,9.231,4.574,3.105,5.805,3.227,6.375,4.153,4.41,4.798,4.838,3.61  /)
    ZTCONSTSTDNIGHT4(1:NCL)=(/ 3.865,4.129,21.824,4.067,69.255,3.922,4.215,10.231,4.279,14.892,6.752,30.25,4.253,4.115,15.248,4.043,17.973,4.151,6.588,40.575,15.923,6.122,4.245,4.114,4.163,4.135,20.778  /)
    ZTCONSTAVGNIGHT5(1:NCL)=(/ 2.669,5.795,4.487,5.643,3.673,4.566,6.701,3.488,5.646,4.346,10.05,9.482,14.929,4.117,3.737,11.956,5.636,8.108,4.877,11.276,3.055,2.747,6.209,4.525,4.492,4.719,5.905  /)
    ZTCONSTSTDNIGHT5(1:NCL)=(/ 17.5,2.635,2.738,2.612,26.686,2.689,13.546,6.808,2.573,6.697,10.153,11.226,54.25,3.019,3.37,25.602,2.635,11.932,2.815,13.72,19.545,11.819,2.737,7.693,2.661,35.406,2.579  /)
    ZTCONSTAVGNIGHT5M(1:NCL)=(/ 2.112,4.58,2.306,4.507,7.92,1.467,4.258,4.919,4.136,9.438,4.776,10.675,4.825,3.929,0.386,10.917,9.521,8.647,4.677,4.473,1.437,6.284,5.644,7.249,3.052,3.652,2.589  /)
    ZTCONSTSTDNIGHT5M(1:NCL)=(/ 2.976,1.554,12.35,1.025,3.07,1.249,11.871,0.865,1.761,5.63,1.311,19.334,25.419,1.17,1.529,7.483,12.356,3.808,1.069,22.031,15.599,1.669,1.158,2.394,1.174,2.699,1.282  /)

! 0.3 Define constants for multi-linear regression of temp 1st layer (only ! dsn<0.20):
ZTMLRIDAY2(1:NCL)=(/2273.776,231.424,-1140.277,1077.753,-2605.152,4792.979,3877.053,393.788,&
                  &-9534.744,-11849.607,748.751,213.009,621.934,-296.89,2569.689,-870.507,&
                  &642.106,775.133,18059.695,218.558,-31549.314,-2667.837,1528.997,&
                  &-1178.686,176.003,2903.15,-2404.29 /)
ZTMLRADAY2(1:NCL)=(/-21.912,-22.981,-18.611,-3.775,19.878,-23.781,-28.105,-6.132,-39.31,&
                  &9.863,-42.609,-27.523,-7.325,-4.131,-21.878,-15.574,7.981,-6.876,5.878,&
                  &-0.976,7.58,-53.716,-14.317,-3.881,0.0,-19.631,-22.49 /)
ZTMLRBDAY2(1:NCL)=(/0.003,-0.013,0.002,0.003,0.001,0.0,-0.008,-0.001,0.003,0.001,-0.003,&
                  &0.003,0.011,-0.002,0.002,0.005,0.008,0.0,0.001,-0.0,0.001,-0.004,0.001,&
                  &0.002,-0.0,-0.004,0.001 /)
ZTMLRCDAY2(1:NCL)=(/2.63,21.656,23.678,-3.807,-0.0,-15.106,0.118,3.81,105.791,79.674,37.359,&
                  &21.693,1.89,7.147,0.852,18.914,-10.844,1.569,-136.832,0.277,225.808,73.195,&
                  &0.27,12.696,0.0,-0.986,41.659 /)
ZTMLRDDAY2(1:NCL)=(/0.046,0.043,0.043,0.007,-0.035,0.052,0.052,0.011,0.081,-0.018,0.079,0.061,&
                  &0.016,0.007,0.045,0.038,-0.019,0.013,-0.014,0.003,-0.015,0.101,0.034,&
                  &0.008,-0.0,0.036,0.041 /)
ZTMLREDAY2(1:NCL)=(/-0.002,-0.037,-0.041,0.01,0.001,0.031,0.003,-0.004,-0.192,-0.148,-0.066,&
                  &-0.037,-0.0,-0.01,0.002,-0.033,0.023,-0.0,0.254,0.0,-0.413,-0.133,0.002,&
                  &-0.021,0.001,0.004,-0.076/)

ZTMLRINIGHT2(1:NCL)=(/25.363,-4915.228,-272.045,-10385.244,-749.825,934.962,-1360.875,246.995,&
                    &-1131.312,1656.874,1945.898,-4666.717,3152.294,1125.078,15549.324,&
                    &-8619.04,-27203.797,769.829,-540.542,521.518,6036.785,127.023,&
                    &6869.104,-1328.559,885.202,-368.013,7971.363 /)
ZTMLRANIGHT2(1:NCL)=(/-21.853,-3.088,-43.653,-29.03,-16.206,1.009,4.797,6.269,-13.664,&
                    &-2.686,24.144,-11.438,-3.761,9.435,-10.072,5.866,-22.161,1.538,&
                    &-27.429,5.321,-30.082,-8.239,-7.267,-32.86,2.993,-41.127,-7.185 /)
ZTMLRBNIGHT2(1:NCL)=(/0.0,0.008,-0.007,0.002,-0.0,0.003,0.009,0.003,0.001,0.004,-0.019,&
                    &0.012,-0.001,0.002,-0.001,0.001,0.0,0.002,-0.002,-0.0,0.006,0.002,&
                    &-0.0,-0.008,-0.001,-0.004,0.0 /)
ZTMLRCNIGHT2(1:NCL)=(/21.915,40.216,45.865,103.003,22.934,-7.261,7.449,-6.425,20.528,&
                    &-9.486,-36.946,44.819,-19.352,-15.056,-103.585,59.556,220.627,-6.459,&
                    &31.715,-7.285,-18.251,7.928,-44.393,43.387,-8.759,43.923,-51.684 /)
ZTMLRDNIGHT2(1:NCL)=(/0.042,0.006,0.082,0.06,0.03,-0.003,-0.013,-0.014,0.029,0.007,-0.05,&
                    &0.026,0.007,-0.023,0.019,-0.01,0.045,-0.004,0.052,-0.013,0.064,0.015,&
                    &0.017,0.061,-0.007,0.078,0.014 /)
ZTMLRENIGHT2(1:NCL)=(/-0.039,-0.073,-0.083,-0.187,-0.041,0.017,-0.012,0.014,-0.035,0.02,0.074,&
                    &-0.081,0.039,0.03,0.192,-0.11,-0.403,0.016,-0.057,0.016,0.037,-0.012,&
                    &0.085,-0.078,0.02,-0.079,0.098/)

ZTMLRIDAY3(1:NCL)=(/-873.735,16851.732,1769.369,-1208.98,595.977,892.861,9584.742,&
                  &515.852,1914.724,-6049.516,-1895.017,-389.405,1705.265,1272.757,&
                  &29055.08,947.23,1665.529,5029.054,-333.738,24368.355,-1959.943,&
                  &1115.958,55.659,-13399.817,857.495,-1397.23,5746.104 /)
ZTMLRADAY3(1:NCL)=(/10.619,-15.266,-10.879,32.929,-25.726,14.589,-17.776,13.11,-27.152,&
                  &-42.891,14.106,-13.758,33.093,-5.929,19.702,40.24,-0.844,-20.965,&
                  &14.592,18.751,-18.933,-42.283,-39.847,7.549,21.17,-4.047,-9.189 /)
ZTMLRBDAY3(1:NCL)=(/0.0,0.0,0.0,0.001,-0.005,-0.008,-0.006,0.002,-0.004,-0.004,0.0,-0.001,&
                  &-0.003,0.004,0.002,-0.008,0.001,0.003,-0.01,-0.001,0.001,-0.008,&
                  &-0.003,0.001,-0.013,-0.001,0.0 /)
ZTMLRCDAY3(1:NCL)=(/-3.02,-109.68,-2.124,-22.655,21.564,-19.85,-53.218,-15.717,12.086,&
                  &84.146,-0.053,17.464,-44.619,-3.168,-232.189,-45.785,-11.335,-17.445,&
                  &-10.411,-196.785,32.479,33.39,39.859,90.671,-25.283,16.025,-33.291 /)
ZTMLRDDAY3(1:NCL)=(/-0.019,0.029,0.021,-0.059,0.049,-0.03,0.034,-0.027,0.052,0.087,-0.025,&
                  &0.026,-0.06,0.012,-0.039,-0.077,0.002,0.041,-0.033,-0.036,0.038,0.081,&
                  &0.075,-0.014,-0.046,0.008,0.018 /)
ZTMLREDAY3(1:NCL)=(/0.007,0.206,0.006,0.041,-0.038,0.041,0.101,0.033,-0.019,-0.153,0.002,&
                  &-0.031,0.083,0.008,0.431,0.088,0.024,0.036,0.025,0.365,-0.058,-0.06,&
                  &-0.072,-0.162,0.052,-0.029,0.064 /)

ZTMLRINIGHT3(1:27)=(/-233.194,1311.565,554.687,-170.099,685.766,5874.629,7978.278,&
                    &-8229.753,626.97,-163.629,1671.225,1383.694,22954.838,-6401.198,&
                    &1016.932,15437.012,465.208,5006.659,231532.438,-4.443,1691.159,768.856,&
                    &-1382.435,-7084.328,7818.812,1345.252,-3449.839 /)
ZTMLRANIGHT3(1:NCL)=(/-22.918,17.157,-22.351,-3.505,15.358,-58.272,-53.205,-3.579,31.338,&
                     &-5.707,15.969,23.453,28.833,-54.761,23.377,-32.093,-3.417,17.97,&
                     &11.031,-13.668,10.341,7.027,11.562,4.886,-1.373,19.383,32.732 /)
ZTMLRBNIGHT3(1:NCL)=(/-0.008,-0.002,-0.006,-0.002,-0.007,0.003,-0.008,-0.003,-0.002,&
                     &-0.007,0.009,-0.004,-0.004,-0.003,-0.006,-0.003,-0.007,-0.001,&
                     &0.017,-0.006,0.022,0.007,0.0,-0.006,-0.009,0.0,0.007 /)
ZTMLRCNIGHT3(1:NCL)=(/25.269,-25.646,16.266,5.886,-19.287,11.638,-8.116,66.062,-32.666,&
                    &7.912,-27.288,-32.353,-194.61,100.024,-29.553,-83.072,0.738,-53.752,&
                    &-1712.245,14.486,-22.205,-11.763,0.009,49.625,-55.516,-28.058,-2.126 /)
ZTMLRDNIGHT3(1:NCL)=(/0.043,-0.034,0.047,0.007,-0.029,0.115,0.104,0.006,-0.064,0.011,-0.032,&
                    &-0.045,-0.059,0.106,-0.045,0.061,0.007,-0.032,-0.021,0.026,-0.021,&
                    &-0.014,-0.021,-0.011,0.002,-0.038,-0.072 /)
ZTMLRENIGHT3(1:NCL)=(/-0.045,0.051,-0.028,-0.01,0.038,-0.019,0.018,-0.121,0.063,-0.013,&
                    &0.055,0.063,0.361,-0.182,0.058,0.158,0.0,0.099,3.15,-0.025,0.045,0.024,&
                    &0.001,-0.09,0.105,0.055,0.008 /)

! 0.4 Define centroids for Density: no distinction day and night...
! A: soT-skt B: soT-snT ; C: rsn/rsnmin
    ZRCENTRA2(1:NCL)=(/ -3541.865,244.7844,-344.0662,426.6829,36.5035,-219.7778,446.1849,69.5413,367.3376,-27.8976,257.3477,-2717.5598,284.9772,&
            & -486.9474,305.8375,-2483.176,-314.6438,-128.4454,-152.9171,438.5479,135.3359,292.2806,37.3985,507.4431,-4451.877,-2862.4211,-132.4309 /)
    ZRCENTRB2(1:NCL)=(/ -102.3517,-587.8513,84.9442,-324.6709,-37.6875,242.4561,-460.9119,-343.4314,45.5019,-306.2778,-380.4187,-70.9393,-158.3012,&
            & -106.4662,-39.5139,-91.8077,-53.4339,-131.6157,39.6053,-461.2061,-75.1291,-213.1413,-500.9586,-193.469,-102.2597,-91.0352,17.6585 /)
    ZRCENTRC2(1:NCL)=(/ -125.8568,424.709,-70.0107,289.3643,-47.6114,-153.1688,359.9082,254.8379,34.638,218.8475,281.4551,-88.5215,162.3987,&
            & 53.5059,-50.0578,-113.7225,35.7485,89.3168,-21.8818,379.2143,-99.3011,194.2758,329.7064,220.6259,-121.7691,-111.6626,-15.906 /)
    ZRCENTRA3(1:NCL)=(/ -135.8665,-17.0828,429.4104,402.9216,-207.8142,-144.7099,-1.596,-111.7997,-107.7011,403.1025,-508.7413,-5289.2378,115.0201,&
            & -507.9801,-5578.6562,-883.5401,-3891.072,-590.0021,153.5239,441.6976,-7720.7236,418.2652,473.0286,-423.3754,47.3167,18.4711,-3023.9482 /)
    ZRCENTRB3(1:NCL)=(/ -721.6349,-582.6187,-145.8941,-190.7771,-59.8339,-451.8439,-512.7739,-257.787,-506.0475,-263.0001,1051.515,-227.4999,-534.728,&
            & 15.5551,-67.8342,-31.5683,-156.2351,-152.8403,-96.3175,-210.2662,-204.5721,-258.3867,-539.0995,-462.1386,-446.4355,-600.0053,-159.447 /)
    ZRCENTRC3(1:NCL)=(/ 446.4028,422.3326,188.6381,227.2264,-78.5153,304.9808,365.7172,107.7415,383.5846,294.4718,-587.9872,-277.8058,430.1257,&
            & -85.229,-87.5357,-43.5049,-196.9755,63.5231,-124.7519,243.4263,-251.0924,290.7506,444.5835,250.0695,338.4489,430.5938,-200.5097 /)
    ZRCENTRA4(1:NCL)=(/ 576.3396,-1919.6548,-263.5831,169.8809,1193.1631,-4648.7017,262.7046,-371.5281,684.8535,-572.1608,867.971,-695.9773,-132.8173,&
            & 210.2056,626.4426,428.1664,-1059.6406,303.5695,-573.3696,-89.3266,-531.5976,-1123.296,663.4559,-143.3086,-1854.3871,-173.9566,766.6892 /)
    ZRCENTRB4(1:NCL)=(/ -1298.3723,-226.0219,-697.0114,-145.742,-836.9321,-262.4575,-1094.834,-593.2739,-511.8727,-47.2088,-827.3694,-544.185,-597.2844,&
            & -43.4083,-1016.2291,53.5328,-279.3172,1093.9867,-459.2903,-709.3953,-41.3843,-64.4666,-38.432,-690.1353,-162.8094,-1053.2731,-1148.8083 /)
    ZRCENTRC4(1:NCL)=(/ 997.8715,-291.5672,459.8056,-194.9653,678.8568,-320.1779,767.1482,372.0334,457.3047,-61.451,685.3994,289.1407,460.6511,&
            & -57.7315,770.9725,79.9255,105.6104,-613.095,283.0966,512.6927,-54.5425,-89.1455,-55.9246,475.4332,-211.9415,744.972,913.0076 /)
    ZRCENTRA5(1:NCL)=(/ 1824.1656,-602.4853,1135.1781,-4194.6084,-646.3918,1503.5759,606.6608,1491.4517,-428.5126,-506.2961,-252.2635,-2836.1387,1466.1106,&
            & -588.9854,-1574.7958,-223.5699,-743.1916,-743.4153,-685.0393,-65.7267,8.6606,-1850.3248,-95.0798,2111.3132,-59.7443,-299.5717,-3411.9668 /)
    ZRCENTRB5(1:NCL)=(/ -4.5079,-1021.686,-2750.4539,-215.6783,-1281.6517,-2953.0793,-2313.7065,-1965.9337,-64.9285,-1129.3881,-1425.0013,-154.5478,-2507.9734,&
            & -1271.6472,-112.3157,-793.4882,-2029.1826,-1273.7791,-1000.6006,-1621.2335,-69.3525,-138.8389,-1829.412,-1814.5347,-2046.4401,-1588.4868,-138.3192 /)
    ZRCENTRC5(1:NCL)=(/ -6.8308,589.8851,1823.7694,-267.1755,808.7336,1971.2648,1504.1857,1349.041,-89.9666,674.03,921.1858,-209.7638,1683.8069,&
            & 714.3266,-154.6993,403.6999,1289.4042,744.8958,597.5707,1039.3883,-92.6542,-188.6363,1201.9548,1250.0219,1334.8756,1018.6988,-185.958 /)
    ZRCENTRA5M(1:NCL)=(/ -414.7569,-619.006,-1992.8358,-2534.5632,-599.2324,-1316.8163,-1476.0657,-1028.012,-1218.9622,-466.2645,-1226.2495,-4905.8696,-2005.9907,&
            & -2037.4348,-915.9469,-1133.9756,-1675.7104,-1589.6157,-1072.4772,-2841.8767,-440.3358,-1543.9376,-840.2651,-950.2744,-1312.647,-531.6243,-2631.9937 /)
    ZRCENTRB5M(1:NCL)=(/ -96.9542,-998.4323,-1540.2692,-980.502,-1216.0809,-65.4584,-1511.527,-1870.0996,-73.6153,-1623.2203,-2237.2195,-1421.1248,-52.5641,&
            & -1732.2773,-1539.8936,-1498.8727,-91.1985,-209.0567,-1921.7404,-1686.1198,-1940.4797,-156.8293,-119.5914,-1102.2733,-1173.5461,-1075.4229,-48.6927 /)
    ZRCENTRC5M(1:NCL)=(/ -150.2658,441.1414,837.5116,216.2132,605.5161,-103.0905,809.4069,1031.671,-115.6356,923.1359,1156.9117,545.6916,-82.0775,&
            & 780.2647,816.7612,803.4692,-142.6787,-272.9421,934.514,737.6534,1131.9098,-232.7525,-180.9293,523.5137,594.5945,519.3375,-75.8856 /)


    ZRCONSTAVG2(1:NCL)=(/ 0.12,0.125,0.126,0.128,0.128,0.132,0.127,0.129,0.129,0.129,0.13,0.127,0.126,&
            & 0.126,0.126,0.115,0.128,0.126,0.126,0.121,0.135,0.123,0.122,0.126,0.111,0.129,0.13 /)
    ZRCONSTSTD2(1:NCL)=(/ 0.101,0.108,0.111,0.112,0.119,0.113,0.108,0.115,0.098,0.112,0.105,0.113,0.109,&
            & 0.107,0.113,0.079,0.107,0.107,0.108,0.112,0.128,0.102,0.1,0.108,0.076,0.113,0.114 /)
    ZRCONSTAVG3(1:NCL)=(/ 4.819,2.627,1.98,2.261,2.97,3.798,3.016,3.261,3.477,1.782,2.638,1.384,3.703,&
            & 3.201,2.413,2.912,2.038,4.479,2.016,1.755,1.237,1.944,2.517,3.891,2.3,3.436,4.548 /)
    ZRCONSTSTD3(1:NCL)=(/ 15.314,11.046,4.718,8.5,6.524,14.134,11.327,8.089,14.615,4.604,5.091,12.03,16.932,&
            & 12.693,5.413,7.349,13.899,16.785,4.696,4.746,4.101,4.381,4.565,13.617,7.62,13.648,22.636 /)
    ZRCONSTAVG4(1:NCL)=(/ 0.686,13.91,0.229,3.642,-2.578,5.589,-0.58,1.08,-1.53,8.571,-1.7,2.769,-0.295,&
            & 17.677,-1.499,-1.369,0.249,-1.665,0.152,0.241,14.771,0.418,-0.335,1.146,11.898,-0.836,-0.306 /)
    ZRCONSTSTD4(1:NCL)=(/ 23.299,50.893,21.59,28.504,3.722,39.499,17.836,24.857,6.354,46.597,12.28,31.879,20.202,&
            & 62.911,12.411,3.432,22.868,15.655,21.344,23.611,56.495,19.416,20.296,27.389,44.558,16.196,20.341 /)
    ZRCONSTAVG5(1:NCL)=(/ 1.218,0.549,0.974,5.675,0.932,0.357,0.881,1.266,3.259,0.713,0.963,6.004,1.178,&
            & 0.804,4.967,1.698,2.387,0.312,0.832,1.658,3.249,5.924,0.988,1.002,1.523,1.143,6.439 /)
    ZRCONSTSTD5(1:NCL)=(/ 1.404,5.054,2.532,29.32,2.878,2.444,3.049,1.805,9.36,4.627,3.077,25.99,2.111,&
            & 3.934,19.76,6.083,11.732,3.503,3.554,1.684,4.874,23.519,3.313,1.623,1.855,3.654,26.931 /)
    ZRCONSTAVG5M(1:NCL)=(/ 0.14,-1.558,-1.342,3.167,-1.077,-0.955,-1.341,-1.134,-0.844,-1.306,1.824,2.787,-1.38,&
            & 2.919,-1.137,-1.165,-1.225,-0.183,2.836,3.391,-1.42,1.358,0.111,-1.651,-1.17,-1.533,-1.179 /)
    ZRCONSTSTD5M(1:NCL)=(/ 2.218,2.379,1.631,18.713,5.933,2.835,1.823,3.616,2.124,1.951,11.624,16.115,1.773,&
            & 16.727,3.866,2.375,1.313,19.134,15.746,18.02,1.967,7.441,3.298,2.649,1.963,3.221,2.22 /)


    ZRMLRI2(1:NCL)=(/ -3541.865,244.7844,-344.0662,426.6829,36.5035,-219.7778,446.1849,69.5413,367.3376,-27.8976,257.3477,-2717.5598,284.9772,&
            & -486.9474,305.8375,-2483.176,-314.6438,-128.4454,-152.9171,438.5479,135.3359,292.2806,37.3985,507.4431,-4451.877,-2862.4211,-132.4309 /)
    ZRMLRA2(1:NCL)=(/ -102.3517,-587.8513,84.9442,-324.6709,-37.6875,242.4561,-460.9119,-343.4314,45.5019,-306.2778,-380.4187,-70.9393,-158.3012,&
            & -106.4662,-39.5139,-91.8077,-53.4339,-131.6157,39.6053,-461.2061,-75.1291,-213.1413,-500.9586,-193.469,-102.2597,-91.0352,17.6585 /)
    ZRMLRB2(1:NCL)=(/ -0.3176,1.0335,0.7129,1.6837,0.9461,0.7601,0.7895,0.5153,3.0846,1.3131,0.3298,-0.4786,1.5664,&
            & 0.109,0.4356,0.7629,0.1642,0.8584,0.5072,4.0052,0.7006,1.6226,0.8063,2.123,2.0176,0.6233,0.8563 /)
    ZRMLRC2(1:NCL)=(/ 14.2216,-0.1808,1.3231,-1.6113,0.0002,0.5571,-1.0224,0.3674,-2.4079,0.4052,-0.1948,11.1977,-1.2242,&
            & 2.3541,-0.6866,9.676,1.5549,0.7283,0.7062,-2.2464,-0.0003,-1.1566,0.666,-2.3386,15.7641,11.0971,0.5535 /)
    ZRMLRD2(1:NCL)=(/ 0.0016,-0.0004,0.0002,-0.0017,-0.0001,0.0004,0.0001,0.0007,-0.0055,-0.0012,0.001,0.0019,-0.0014,&
            & 0.0014,0.0007,0.0002,0.0016,-0.0001,0.0007,-0.0085,0.0002,-0.0018,-0.0,-0.0027,-0.0011,0.0004,-0.0001 /)
    ZRMLRE2(1:NCL)=(/ -125.8568,424.709,-70.0107,289.3643,-47.6114,-153.1688,359.9082,254.8379,34.638,218.8475,281.4551,-88.5215,162.3987,&
            & 53.5059,-50.0578,-113.7225,35.7485,89.3168,-21.8818,379.2143,-99.3011,194.2758,329.7064,220.6259,-121.7691,-111.6626,-15.906 /)
    ZRMLRI3(1:NCL)=(/ -135.8665,-17.0828,429.4104,402.9216,-207.8142,-144.7099,-1.596,-111.7997,-107.7011,403.1025,-508.7413,-5289.2378,115.0201,&
            & -507.9801,-5578.6562,-883.5401,-3891.072,-590.0021,153.5239,441.6976,-7720.7236,418.2652,473.0286,-423.3754,47.3167,18.4711,-3023.9482 /)
    ZRMLRA3(1:NCL)=(/ -721.6349,-582.6187,-145.8941,-190.7771,-59.8339,-451.8439,-512.7739,-257.787,-506.0475,-263.0001,1051.515,-227.4999,-534.728,&
            & 15.5551,-67.8342,-31.5683,-156.2351,-152.8403,-96.3175,-210.2662,-204.5721,-258.3867,-539.0995,-462.1386,-446.4355,-600.0053,-159.447 /)
    ZRMLRB3(1:NCL)=(/ -0.486,0.5566,4.1942,2.6247,0.2497,0.5331,0.6235,-0.7123,0.8447,3.2066,5.8237,-1.074,1.4969,&
            & -0.7893,0.3735,0.3227,-1.1986,0.2935,0.6345,3.6235,-3.1071,3.0857,5.8253,-0.2094,0.773,1.0088,0.0695 /)
    ZRMLRC3(1:NCL)=(/ 2.1963,0.9702,-2.8915,-2.1389,1.3874,1.3616,0.8246,1.6283,1.0291,-2.3191,-1.6234,21.6286,-0.0761,&
            & 2.7677,21.0146,3.7256,16.3228,2.7988,-0.0002,-2.6451,31.8912,-2.3259,-2.9774,2.8578,0.4271,0.7221,12.3174 /)
    ZRMLRD3(1:NCL)=(/ 0.0023,0.0004,-0.0081,-0.0041,0.0008,0.0002,0.0002,0.0027,0.0,-0.0054,-0.0125,0.0025,-0.0015,&
            & 0.0028,0.0006,0.0008,0.0027,0.0006,0.0002,-0.0067,0.005,-0.0052,-0.0134,0.0016,0.0001,-0.0009,0.0011 /)
    ZRMLRE3(1:NCL)=(/ 446.4028,422.3326,188.6381,227.2264,-78.5153,304.9808,365.7172,107.7415,383.5846,294.4718,-587.9872,-277.8058,430.1257,&
            & -85.229,-87.5357,-43.5049,-196.9755,63.5231,-124.7519,243.4263,-251.0924,290.7506,444.5835,250.0695,338.4489,430.5938,-200.5097 /)
    ZRMLRI4(1:NCL)=(/ 576.3396,-1919.6548,-263.5831,169.8809,1193.1631,-4648.7017,262.7046,-371.5281,684.8535,-572.1608,867.971,-695.9773,-132.8173,&
            & 210.2056,626.4426,428.1664,-1059.6406,303.5695,-573.3696,-89.3266,-531.5976,-1123.296,663.4559,-143.3086,-1854.3871,-173.9566,766.6892 /)
    ZRMLRA4(1:NCL)=(/ -1298.3723,-226.0219,-697.0114,-145.742,-836.9321,-262.4575,-1094.834,-593.2739,-511.8727,-47.2088,-827.3694,-544.185,-597.2844,&
            & -43.4083,-1016.2291,53.5328,-279.3172,1093.9867,-459.2903,-709.3953,-41.3843,-64.4666,-38.432,-690.1353,-162.8094,-1053.2731,-1148.8083 /)
    ZRMLRB4(1:NCL)=(/ 2.7242,-2.701,0.375,0.995,2.2678,-1.556,0.8805,0.1763,2.2016,-0.7607,2.6228,-0.7007,1.2143,&
            & -0.416,1.3524,1.9215,-0.3578,2.6048,0.1285,1.2984,-0.4526,-0.0033,0.5199,0.3166,-0.8742,0.737,2.7089 /)
    ZRMLRC4(1:NCL)=(/ -1.3924,10.3721,2.1816,-0.0005,-4.2505,19.6819,0.4233,2.5894,-2.7521,3.125,-3.1591,4.1866,1.053,&
            & 0.0772,-1.4036,-2.5888,5.0473,-3.7766,3.1704,1.0984,2.8336,4.8803,-2.2464,1.6863,8.5633,2.0769,-2.3699 /)
    ZRMLRD4(1:NCL)=(/ -0.0043,0.0047,0.0006,-0.0005,-0.0031,0.0032,-0.0001,0.0008,-0.0029,0.0026,-0.004,0.0026,-0.0009,&
            & 0.002,-0.0009,-0.0013,0.0018,-0.0038,0.0009,-0.0014,0.0019,0.0015,0.0009,0.0009,0.0024,0.0001,-0.0042 /)
    ZRMLRE4(1:NCL)=(/ 997.8715,-291.5672,459.8056,-194.9653,678.8568,-320.1779,767.1482,372.0334,457.3047,-61.451,685.3994,289.1407,460.6511,&
            & -57.7315,770.9725,79.9255,105.6104,-613.095,283.0966,512.6927,-54.5425,-89.1455,-55.9246,475.4332,-211.9415,744.972,913.0076 /)
    ZRMLRI5(1:NCL)=(/ 1824.1656,-602.4853,1135.1781,-4194.6084,-646.3918,1503.5759,606.6608,1491.4517,-428.5126,-506.2961,-252.2635,-2836.1387,1466.1106,&
            & -588.9854,-1574.7958,-223.5699,-743.1916,-743.4153,-685.0393,-65.7267,8.6606,-1850.3248,-95.0798,2111.3132,-59.7443,-299.5717,-3411.9668 /)
    ZRMLRA5(1:NCL)=(/ -4.5079,-1021.686,-2750.4539,-215.6783,-1281.6517,-2953.0793,-2313.7065,-1965.9337,-64.9285,-1129.3881,-1425.0013,-154.5478,-2507.9734,&
            & -1271.6472,-112.3157,-793.4882,-2029.1826,-1273.7791,-1000.6006,-1621.2335,-69.3525,-138.8389,-1829.412,-1814.5347,-2046.4401,-1588.4868,-138.3192 /)
    ZRMLRB5(1:NCL)=(/ -0.1663,-0.0223,2.6378,-0.23,0.8619,2.7671,2.5698,1.7874,-1.3583,0.3275,1.0704,-1.8988,2.3386,&
            & -1.0912,-1.3509,-1.5455,0.0438,0.0389,0.4111,1.687,-1.7107,-1.6032,1.6402,2.2335,1.4764,1.1552,-1.9405 /)
    ZRMLRC5(1:NCL)=(/ -6.7465,4.3331,-1.1129,16.7992,4.3355,-2.2735,0.4202,-3.439,2.9963,3.9399,2.9283,12.7039,-2.6644,&
            & 5.1551,7.4995,3.3386,6.0152,5.2089,4.3586,2.2071,1.5721,8.7876,2.6435,-6.2158,2.8229,3.304,14.729 /)
    ZRMLRD5(1:NCL)=(/ 0.0026,0.0006,-0.0047,0.0018,-0.001,-0.0051,-0.0049,-0.0024,0.0035,-0.0001,-0.0013,0.004,-0.004,&
            & 0.0029,0.0034,0.0038,0.0011,0.0004,-0.0002,-0.0025,0.0042,0.0037,-0.0025,-0.0035,-0.0019,-0.0013,0.0042 /)
    ZRMLRE5(1:NCL)=(/ -6.8308,589.8851,1823.7694,-267.1755,808.7336,1971.2648,1504.1857,1349.041,-89.9666,674.03,921.1858,-209.7638,1683.8069,&
            & 714.3266,-154.6993,403.6999,1289.4042,744.8958,597.5707,1039.3883,-92.6542,-188.6363,1201.9548,1250.0219,1334.8756,1018.6988,-185.958 /)
    ZRMLRI5M(1:NCL)=(/ -414.7569,-619.006,-1992.8358,-2534.5632,-599.2324,-1316.8163,-1476.0657,-1028.012,-1218.9622,-466.2645,-1226.2495,-4905.8696,-2005.9907,&
            & -2037.4348,-915.9469,-1133.9756,-1675.7104,-1589.6157,-1072.4772,-2841.8767,-440.3358,-1543.9376,-840.2651,-950.2744,-1312.647,-531.6243,-2631.9937 /)
    ZRMLRA5M(1:NCL)=(/ -96.9542,-998.4323,-1540.2692,-980.502,-1216.0809,-65.4584,-1511.527,-1870.0996,-73.6153,-1623.2203,-2237.2195,-1421.1248,-52.5641,&
            & -1732.2773,-1539.8936,-1498.8727,-91.1985,-209.0567,-1921.7404,-1686.1198,-1940.4797,-156.8293,-119.5914,-1102.2733,-1173.5461,-1075.4229,-48.6927 /)
    ZRMLRB5M(1:NCL)=(/ -0.3551,-0.5878,1.1112,-3.0124,0.3306,3.511,1.2367,0.4404,2.0579,0.3505,-2.231,-3.2196,1.125,&
            & -2.8588,0.3113,1.1335,6.1145,-1.6367,-2.742,-2.8471,0.6615,-1.4675,-0.5128,-0.7918,1.7464,0.1495,1.3638 /)
    ZRMLRC5M(1:NCL)=(/ 2.9368,4.9979,9.982,13.7965,4.7362,4.2015,8.0027,7.1985,4.5569,4.673,10.0353,23.1413,7.7932,&
            & 12.7489,6.3764,6.7421,4.7105,8.1889,9.3455,15.6609,4.8956,7.9932,4.7174,6.4062,6.663,4.3699,9.9766 /)
    ZRMLRD5M(1:NCL)=(/ 0.0006,0.0012,-0.002,0.0052,-0.0005,-0.0068,-0.0023,-0.0007,-0.0038,-0.0004,0.004,0.0055,-0.002,&
            & 0.0051,-0.0004,-0.0021,-0.0128,0.0037,0.0048,0.005,-0.0011,0.0027,0.0009,0.0017,-0.0033,-0.0002,-0.0025 /)
    ZRMLRE5M(1:NCL)=(/ -150.2658,441.1414,837.5116,216.2132,605.5161,-103.0905,809.4069,1031.671,-115.6356,923.1359,1156.9117,545.6916,-82.0775,&
            & 780.2647,816.7612,803.4692,-142.6787,-272.9421,934.514,737.6534,1131.9098,-232.7525,-180.9293,523.5137,594.5945,519.3375,-75.8856 /)


!*******************************************************************************
! 2. Start snow parametrizations
!    Here we start the snow temperature and density parametrisation:
!    both are simple exponential function. Temperature relaxes to first soil layer
!    at the bottom, while the density relaxes to mean density. 
!    The missing snow mass it is added to the bottom layer simply increasing the
!    density of this layer keeping the depth fixed.
!*******************************************************************************
  DO JL=KIDIA, KFDIA
    IF (PMU0(JL) > ZEPSILON ) THEN
      IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY2
        ZTCENTRB=ZTCENTRBDAY2
        ZTCENTRC=ZTCENTRCDAY2

        ZTCONSTAVG=ZTCONSTAVGDAY2
        ZTCONSTSTD=ZTCONSTSTDDAY2

        ZTMLRI=ZTMLRIDAY2
        ZTMLRA=ZTMLRADAY2
        ZTMLRB=ZTMLRBDAY2
        ZTMLRC=ZTMLRCDAY2
        ZTMLRD=ZTMLRDDAY2
        ZTMLRE=ZTMLREDAY2

      ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY3
        ZTCENTRB=ZTCENTRBDAY3
        ZTCENTRC=ZTCENTRCDAY3

        ZTCONSTAVG=ZTCONSTAVGDAY3
        ZTCONSTSTD=ZTCONSTSTDDAY3
        
        ZTMLRI=ZTMLRIDAY3
        ZTMLRA=ZTMLRADAY3
        ZTMLRB=ZTMLRBDAY3
        ZTMLRC=ZTMLRCDAY3
        ZTMLRD=ZTMLRDDAY3
        ZTMLRE=ZTMLREDAY3

      ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY4
        ZTCENTRB=ZTCENTRBDAY4
        ZTCENTRC=ZTCENTRCDAY4
        
        ZTCONSTAVG=ZTCONSTAVGDAY4
        ZTCONSTSTD=ZTCONSTSTDDAY4

      ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY5
        ZTCENTRB=ZTCENTRBDAY5
        ZTCENTRC=ZTCENTRCDAY5

        ZTCONSTAVG=ZTCONSTAVGDAY5
        ZTCONSTSTD=ZTCONSTSTDDAY5
        
      ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRADAY5M
        ZTCENTRB=ZTCENTRBDAY5M
        ZTCENTRC=ZTCENTRCDAY5M

        ZTCONSTAVG=ZTCONSTAVGDAY5M
        ZTCONSTSTD=ZTCONSTSTDDAY5M
      ENDIF
        
    ELSEIF (PMU0(JL)<=ZEPSILON) THEN
      IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT2
        ZTCENTRB=ZTCENTRBNIGHT2
        ZTCENTRC=ZTCENTRCNIGHT2

        ZTCONSTAVG=ZTCONSTAVGNIGHT2
        ZTCONSTSTD=ZTCONSTSTDNIGHT2

        ZTMLRI=ZTMLRINIGHT2
        ZTMLRA=ZTMLRANIGHT2
        ZTMLRB=ZTMLRBNIGHT2
        ZTMLRC=ZTMLRCNIGHT2
        ZTMLRD=ZTMLRDNIGHT2
        ZTMLRE=ZTMLRENIGHT2
        
      ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT3
        ZTCENTRB=ZTCENTRBNIGHT3
        ZTCENTRC=ZTCENTRCNIGHT3

        ZTCONSTAVG=ZTCONSTAVGNIGHT3
        ZTCONSTSTD=ZTCONSTSTDNIGHT3
        
        ZTMLRI=ZTMLRINIGHT3
        ZTMLRA=ZTMLRANIGHT3
        ZTMLRB=ZTMLRBNIGHT3
        ZTMLRC=ZTMLRCNIGHT3
        ZTMLRD=ZTMLRDNIGHT3
        ZTMLRE=ZTMLRENIGHT3

      ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT4
        ZTCENTRB=ZTCENTRBNIGHT4
        ZTCENTRC=ZTCENTRCNIGHT4

        ZTCONSTAVG=ZTCONSTAVGNIGHT4
        ZTCONSTSTD=ZTCONSTSTDNIGHT4
        
      ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT5
        ZTCENTRB=ZTCENTRBNIGHT5
        ZTCENTRC=ZTCENTRCNIGHT5

        ZTCONSTAVG=ZTCONSTAVGNIGHT5
        ZTCONSTSTD=ZTCONSTSTDNIGHT5
        
      ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
        ZTCENTRA=ZTCENTRANIGHT5M
        ZTCENTRB=ZTCENTRBNIGHT5M
        ZTCENTRC=ZTCENTRCNIGHT5M

        ZTCONSTAVG=ZTCONSTAVGNIGHT5M
        ZTCONSTSTD=ZTCONSTSTDNIGHT5M
      ENDIF
    ENDIF

    PTCONSTAVG(JL,1:NCL)=ZTCONSTAVG(1:NCL)
    PTCONSTSTD(JL,1:NCL)=ZTCONSTSTD(1:NCL)

! Assign density depending on snow depth::
    IF ( ZDSNTOT(JL) < 0.15_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA2
      ZRCENTRB=ZRCENTRB2
      ZRCENTRC=ZRCENTRC2

      ZRCONSTAVG=ZRCONSTAVG2
      ZRCONSTSTD=ZRCONSTSTD2

      ZRMLRI=ZRMLRI2
      ZRMLRA=ZRMLRA2
      ZRMLRB=ZRMLRB2
      ZRMLRC=ZRMLRC2
      ZRMLRD=ZRMLRD2
      ZRMLRE=ZRMLRE2
      
      KLEVMID(JL)=2_JPIM
    ELSEIF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA3
      ZRCENTRB=ZRCENTRB3
      ZRCENTRC=ZRCENTRC3

      ZRCONSTAVG=ZRCONSTAVG3
      ZRCONSTSTD=ZRCONSTSTD3
      
      ZRMLRI=ZRMLRI3
      ZRMLRA=ZRMLRA3
      ZRMLRB=ZRMLRB3
      ZRMLRC=ZRMLRC3
      ZRMLRD=ZRMLRD3
      ZRMLRE=ZRMLRE3

      KLEVMID(JL)=2_JPIM

    ELSEIF ( ZDSNTOT(JL) < 0.25_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA4
      ZRCENTRB=ZRCENTRB4
      ZRCENTRC=ZRCENTRC4

      ZRCONSTAVG=ZRCONSTAVG4
      ZRCONSTSTD=ZRCONSTSTD4
      
      ZRMLRI=ZRMLRI4
      ZRMLRA=ZRMLRA4
      ZRMLRB=ZRMLRB4
      ZRMLRC=ZRMLRC4
      ZRMLRD=ZRMLRD4
      ZRMLRE=ZRMLRE4
 
      KLEVMID(JL)=3_JPIM

    ELSEIF ( ZDSNTOT(JL) < 0.50_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA5
      ZRCENTRB=ZRCENTRB5
      ZRCENTRC=ZRCENTRC5

      ZRCONSTAVG=ZRCONSTAVG5
      ZRCONSTSTD=ZRCONSTSTD5
      
      ZRMLRI=ZRMLRI5
      ZRMLRA=ZRMLRA5
      ZRMLRB=ZRMLRB5
      ZRMLRC=ZRMLRC5
      ZRMLRD=ZRMLRD5
      ZRMLRE=ZRMLRE5

      IF (PSDOR(JL)<ZSDORTHR)THEN
        KLEVMID(JL)=3_JPIM
      ELSEIF (KLEVSNA(JL)<KLEVSN .AND. KLEVSNA(JL)>2)THEN
        KLEVMID(JL)=KLEVSNA(JL)-1
      ELSEIF (KLEVSNA(JL)==2)THEN
        KLEVMID(JL)=KLEVSNA(JL)
      ELSE
        KLEVMID(JL)=3_JPIM
      ENDIF

    ELSEIF ( ZDSNTOT(JL) >= 0.50_JPRB ) THEN 
      ZRCENTRA=ZRCENTRA5M
      ZRCENTRB=ZRCENTRB5M
      ZRCENTRC=ZRCENTRC5M

      ZRCONSTAVG=ZRCONSTAVG5M
      ZRCONSTSTD=ZRCONSTSTD5M
      
      ZRMLRI=ZRMLRI5M
      ZRMLRA=ZRMLRA5M
      ZRMLRB=ZRMLRB5M
      ZRMLRC=ZRMLRC5M
      ZRMLRD=ZRMLRD5M
      ZRMLRE=ZRMLRE5M

      KLEVMID(JL)=MAX(KLEVSNA(JL)-1,1)
     
    ENDIF
! Final assignment:
!  RCENTRA(JL)=ZRCENTRA
!  RCENTRB(JL)=ZRCENTRB
!  RCENTRC(JL)=ZRCENTRC

   PRCONSTAVG(JL,1:NCL)=ZRCONSTAVG(1:NCL)
   PRCONSTSTD(JL,1:NCL)=ZRCONSTSTD(1:NCL)
   
!  RMLRI(JL)=ZRMLRI
!  RMLRA(JL)=ZRMLRA
!  RMLRB(JL)=ZRMLRB
!  RMLRC(JL)=ZRMLRC
!  RMLRD(JL)=ZRMLRD
!  RMLRE(JL)=ZRMLRE

!*****************************************
! 2.1 Find closest cluster centre using euclidean metric:
!     Features are common to temp and dens
   ZFEATA=PTSOIL(JL)-PTSKIN(JL)
   ZFEATB=PTSOIL(JL)-PTSN(JL)
   ZFEATC=PRSN(JL)/RHOMINSND

! 2.1.1 Temperature
   DO KCL=1,NCL
     ZTDIST(KCL)=SQRT( (ZFEATA-ZTCENTRA(KCL))**2_JPRB + (ZFEATB-ZTCENTRB(KCL))**2_JPRB + (ZFEATC-ZTCENTRC(KCL))**2_JPRB )
   ENDDO
   PTMINCL(JL)=MINLOC(ZTDIST,DIM=1)

! 2.1.2 Density:
   DO KCL=1,NCL
     ZRDIST(KCL)=SQRT( (ZFEATA-ZRCENTRA(KCL))**2_JPRB + (ZFEATB-ZRCENTRB(KCL))**2_JPRB + (ZFEATC-ZRCENTRC(KCL))**2_JPRB )
   ENDDO
   PRMINCL(JL)=MINLOC(ZRDIST,DIM=1)

!*****************************************
! 2.1 Initialize warm start (WS) variables
    PTSNTOP(JL)=PTSKIN(JL)
    IF ( PSSN(JL) < ZSNPERT ) THEN
        PTSNWS(JL, 1)          = PTSKIN(JL)
        PTSNWS(JL, 2:KLEVSN-1) = PTSN(JL)
        PTSNWS(JL, KLEVSN)     = PTSOIL(JL)
        PTSNMIDDLE(JL)         = PTSN(JL)

        PRSNWS(JL, 1:KLEVSN) = PRSN(JL)
        PRSNMAX(JL)          = RHOMAXSN_NEW

        PSSNWS(JL, 1)        = PSSN(JL)
        PWSNWS(JL, 1)        = PWSN(JL)
        PSSNWS(JL, 2:KLEVSN) = 0._JPRB
        PWSNWS(JL, 2:KLEVSN) = 0._JPRB

        IF ( ZDSNTOT(JL) < 0.20_JPRB ) THEN
          PTSNBOTTOM(JL)=PTSOIL(JL)
          ZSADEPTH(JL)=0._JPRB
        ELSE
          ZSADEPTH(JL)=0.5_JPRB*RDAT(1)
          IF (PTSNTOP(JL)>=PTSN(JL)) THEN
            PTSNBOTTOM=PTSN(JL)
           !IF (PSSN(JL) < ZSNPERT) THEN
              ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL)-1)
           !ELSE
           !  ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL))
           !ENDIF
          ELSE
            PTSNBOTTOM=PTSOIL(JL)
            ZACTDEPTH(JL)=ZDSNTOT(JL)
          ENDIF
        ENDIF
!----- Density
      ZP1=ZRMLRI(PRMINCL(JL))+ZRMLRA(PRMINCL(JL))*PASN(JL)
      ZP2=                    ZRMLRB(PRMINCL(JL))*PRSN(JL)
      ZP3=                    ZRMLRC(PRMINCL(JL))*PTSN(JL)
      ZP4=                    ZRMLRD(PRMINCL(JL))*(PRSN(JL))**2._JPRB
      ZP5=                    ZRMLRE(PRMINCL(JL))*(PASN(JL))**2._JPRB

      PRSNTOP(JL)=ZP1+ZP2+ZP3+ZP4+ZP5
      !!!!!
      ! PTSNBOTTOM(JL)=PTSOIL(JL)
      ! PTSNMIDDLE(JL)=0.5*(PTSN(JL)+PTSOIL(JL))
      ! PTSNTOP(JL)   =PTSN(JL)
      !!!!!

    ELSE ! Glacier ini
        KLEVSNA(JL)=KLEVSN !KSNACC+1

        PTSNWS(JL, 1)        = PTSKIN(JL)
        PTSNWS(JL, 2:KLEVSN) = PTSN(JL)
       !PTSNWS(JL, KLEVSN+1) = PTSN(JL)

        PTSNMIDDLE(JL)=PTSN(JL)
        PTSNBOTTOM(JL)=PTSN(JL)

        PRSNMAX(JL)              = 300._JPRB
        PRSNWS(JL,1:KLEVSN)      = PRSNMAX(JL)
       !PRSNWS(JL,KSNACC:KLEVSN) = PRSNMAX(JL)
        PSSNWS(JL,1:KLEVSN)      = PRSNMAX(JL)*ZDSNR(JL,1:KLEVSN)
        PWSNWS(JL,1:KLEVSN)      = 0._JPRB
!----- Active depth
      IF (PTSNTOP(JL)>=PTSN(JL)) THEN
        ZACTDEPTH(JL)=ZSNDEPTH(JL, KLEVSNA(JL))
      ELSE
        ZACTDEPTH(JL)=ZDSNTOT(JL)
      ENDIF
!----- Density
      PRSNTOP(JL)=300._JPRB
    ENDIF
  ENDDO 
END ASSOCIATE

!    -----------------------------------------------------------------
IF (LHOOK) CALL DR_HOOK('SURFWS_INIT_MLOFF_MOD:SURFWS_INIT_MLOFF',1,ZHOOK_HANDLE)

END SUBROUTINE SURFWS_INIT_MLOFF
END MODULE SURFWS_INIT_MLOFF_MOD


